//-------------------------------------------------------------------------------------------//
// Copyright(c) Wang Yong
// Author     : Wang Yong
// Date       : 2021-05-30
// Description: Refine the tetrahedron mesh by split the marked element. It is designed for
// FEM adaptive mesh function.
//-------------------------------------------------------------------------------------------//
#include <cmath>
#include <algorithm>
#include "RefineTetr.h"

constexpr Real EDGE_COMPARE = 1.2;

// find all common edge
void RefineTetr::fn_FindTetrEdges(std::vector<Node> &vR, std::vector<Tetr> &vElem)
{
    constexpr int LocNode[6][2] = { {0,1}, {1,2}, {0,2}, {0,3}, {1,3}, {2,3} };

    if (vElem.empty()) return;

    int iNode, m, n, j, k;
    bool bFind;

    m_vTetrEdge.resize(vElem.size(), std::vector<int>(6));
    auto vvNode2Edge = std::vector< std::vector<int>>(vR.size());
    // set first tetrahedron
    std::vector< std::vector<int>> vvLib;
    auto vEdge = std::vector<int>(2);
    for (auto i = 0; i < 6; ++i) {
        m_vTetrEdge[0][i] = i;
        for (j = 0; j < 2; ++j) {
            iNode = vElem[0].node[LocNode[i][j]];
            vEdge[j] = iNode;
            vvNode2Edge[iNode].push_back(i);
        }
        std::sort(vEdge.begin(), vEdge.end());
        vvLib.push_back(vEdge);
    }
    // check the other tetrahedron
    for (auto i = 1; i < vElem.size(); ++i) {
        // check each edge
        for (n = 0; n < 6; ++n) {
            vEdge[0] = vElem[i].node[LocNode[n][0]];
            vEdge[1] = vElem[i].node[LocNode[n][1]];
            std::sort(vEdge.begin(), vEdge.end());
            bFind = false;
            for (j = 0; j < 2; ++j) {
                iNode = vEdge[j];
                for (k = 0; k < vvNode2Edge[iNode].size(); ++k) {
                    m = vvNode2Edge[iNode][k];
                    if (vvLib[m][0] == vEdge[0] && vvLib[m][1] == vEdge[1]) {
                        bFind = true;
                        m_vTetrEdge[i][n] = m;
                        break;
                    }
                }
                if (bFind) break;
            }
            if (!bFind) {
                m = int(vvLib.size());
                m_vTetrEdge[i][n] = m;
                vvNode2Edge[vEdge[0]].push_back(m);
                vvNode2Edge[vEdge[1]].push_back(m);
                vvLib.push_back(vEdge);
            }
        }
    }
    // get each edge length
    Node r;
    Real dist;
    m_vEdgeLen.reserve(vvLib.size());
    for (auto& v : vvLib) {
        r.x = vR[v[0]].x - vR[v[1]].x;
        r.y = vR[v[0]].y - vR[v[1]].y;
        r.z = vR[v[0]].z - vR[v[1]].z;
        dist = sqrt(r.x * r.x + r.y * r.y + r.z * r.z);
        m_vEdgeLen.push_back(dist);
    }
}

// mark the edges need to be splitted
void RefineTetr::fn_MarkSplitEdge(std::vector<Tetr> &vElem, std::vector<bool> &vRefine)
{
    if (vElem.empty()) return;

    Real dMinEdge, dMaxEdge;
    int iMaxEdge, j;

    m_vEdgeFlag.resize(m_vEdgeLen.size(), false);
    // check all element edge
    for (auto i = 0; i < vElem.size(); ++i) {
        if (vRefine[i]) {
            iMaxEdge = 0;
            dMinEdge = dMaxEdge = m_vEdgeLen[m_vTetrEdge[i][0]];
            for (j = 1; j < 6; ++j) {
                if (dMinEdge > m_vEdgeLen[m_vTetrEdge[i][j]]) {
                    dMinEdge = m_vEdgeLen[m_vTetrEdge[i][j]];
                }
                if (dMaxEdge < m_vEdgeLen[m_vTetrEdge[i][j]]) {
                    dMaxEdge = m_vEdgeLen[m_vTetrEdge[i][j]];
                    iMaxEdge = j;
                }
            }
            dMinEdge *= EDGE_COMPARE;
            dMinEdge = std::min(dMinEdge, dMaxEdge * 0.6);

            m_vEdgeFlag[m_vTetrEdge[i][iMaxEdge]] = true;
            for (j = 0; j < 6; ++j) {
                if (m_vEdgeLen[m_vTetrEdge[i][j]] > dMinEdge) {
                    m_vEdgeFlag[m_vTetrEdge[i][j]] = true;
                }
            }
        }
    }
}

// adapt the edge to improve the mesh quality
void RefineTetr::fn_AdaptTetrEdge(std::vector<Tetr> &vElem)
{
    Real dMinEdge;
    int edge[6], iEdge, m, j;

    for (auto i = 0; i < vElem.size(); ++i) {
        m = 0;
        for (j = 0; j < 6; ++j) {
            if (m_vEdgeFlag[m_vTetrEdge[i][j]]) {
                edge[m] = m_vTetrEdge[i][j];
                ++m;
            }
        }
        if (m > 0) {
            dMinEdge = m_vEdgeLen[edge[0]];
            for (j = 1; j < m; ++j) {
                if (dMinEdge > m_vEdgeLen[edge[j]]) {
                    dMinEdge = m_vEdgeLen[edge[j]];
                }
            }
            for (j = 0; j < 6; ++j) {
                iEdge = m_vTetrEdge[i][j];
                if (m_vEdgeLen[iEdge] > dMinEdge) {
                    m_vEdgeFlag[iEdge] = true;
                }
            }
        }
    }
}

// check and update the split edge number
void RefineTetr::fn_CheckTetrType(std::vector<Tetr> &vElem)
{
    if (vElem.empty()) return;

    int e[6], j, m;

    // set the type of each element
    m_vTetrType.resize(vElem.size(), 0);
    for (auto i = 0; i < vElem.size(); ++i) {
        m = 0;
        for (j = 0; j < 6; ++j) {
            if (m_vEdgeFlag[m_vTetrEdge[i][j]]) {
                e[m] = j;
                ++m;
            }
        }
        switch (m) {
        case(1):
            m_vTetrType[i] = Edge1Type1;
            break;
        case(2):
            if (e[0] == 0 && e[1] == 5) {
                m_vTetrType[i] = Edge2Type2;
            } else if (e[0] == 1 && e[1] == 3) {
                m_vTetrType[i] = Edge2Type2;
            } else if (e[0] == 2 && e[1] == 4) {
                m_vTetrType[i] = Edge2Type2;
            } else {
                m_vTetrType[i] = Edge2Type3;
            }
            break;
        case(3):
            if (e[0] == 0 && e[1] == 1 && e[2] == 2) {
                m_vTetrType[i] = Edge3Type4;
            } else if (e[0] == 1 && e[1] == 4 && e[2] == 5) {
                m_vTetrType[i] = Edge3Type4;
            } else if (e[0] == 2 && e[1] == 3 && e[2] == 5) {
                m_vTetrType[i] = Edge3Type4;
            } else if (e[0] == 0 && e[1] == 3 && e[2] == 4) {
                m_vTetrType[i] = Edge3Type4;
            } else if (e[0] == 0 && e[1] == 2 && e[2] == 3) {
                m_vTetrType[i] = Edge3Type5;
            } else if (e[0] == 0 && e[1] == 1 && e[2] == 4) {
                m_vTetrType[i] = Edge3Type5;
            } else if (e[0] == 1 && e[1] == 2 && e[2] == 5) {
                m_vTetrType[i] = Edge3Type5;
            } else if (e[0] == 3 && e[1] == 4 && e[2] == 5) {
                m_vTetrType[i] = Edge3Type5;
            } else {
                m_vTetrType[i] = Edge3Type6;
            }
            break;
        case(4):
            if (e[0] == 1 && e[1] == 2 && e[2] == 3 && e[3] == 4) {
                m_vTetrType[i] = Edge4Type7;
            } else if (e[0] == 0 && e[1] == 2 && e[2] == 4 && e[3] == 5) {
                m_vTetrType[i] = Edge4Type7;
            } else if (e[0] == 0 && e[1] == 1 && e[2] == 3 && e[3] == 5) {
                m_vTetrType[i] = Edge4Type7;
            } else {
                m_vTetrType[i] = Edge4Type8;
            }
            break;
        case(5):
            m_vTetrType[i] = Edge5Type9;
            break;
        case(6):
            m_vTetrType[i] = AllSixEdge;
            break;
        default:
            m_vTetrType[i] = NOT_SPLIT;
            break;
        }
    }
}

// insert new point on each edge
void RefineTetr::fn_AddEdgesPoint(std::vector<Tetr> &vElem, std::vector<Node> &vR)
{
    constexpr int LocNode[6][2] = { {0,1}, {1,2}, {0,2}, {0,3}, {1,3}, {2,3} };

    if (vElem.empty()) return;

    Node r;
    int iEdge, j;

    // count how many point will be created;
    iEdge = 0;
    for (auto cut : m_vEdgeFlag) {
        if (cut) ++iEdge;
    }
    vR.reserve(vR.size() + iEdge);

    m_vEdgeNode.resize(m_vEdgeLen.size(), -1);
    // check all element edge
    for (auto i = 0; i < vElem.size(); ++i) {
        for (j = 0; j < 6; ++j) {
            iEdge = m_vTetrEdge[i][j];
            if (m_vEdgeFlag[iEdge] && m_vEdgeNode[iEdge] < 0) {
                m_vEdgeNode[iEdge] = int(vR.size());
                r.x = vR[vElem[i].node[LocNode[j][0]]].x + vR[vElem[i].node[LocNode[j][1]]].x;
                r.y = vR[vElem[i].node[LocNode[j][0]]].y + vR[vElem[i].node[LocNode[j][1]]].y;
                r.z = vR[vElem[i].node[LocNode[j][0]]].z + vR[vElem[i].node[LocNode[j][1]]].z;
                r.x *= 0.5;
                r.y *= 0.5;
                r.z *= 0.5;
                vR.push_back(r);
            }
        }
    }
}

// split each element
void RefineTetr::fn_SplitEachTetr(std::vector<Node> &vR, std::vector<Tetr> &vElem)
{
    if (vElem.empty()) return;

    std::vector<Tetr> vNewElem;
    for (auto i = 0; i < vElem.size(); ++i) {
        switch (m_vTetrType[i]) {
        case(Edge1Type1):
            // split only one tetrahedra edge
            fn_CutEdge1Type1(vElem[i], m_vTetrEdge[i], vNewElem);
            break;
        case(Edge2Type2):
            // split two tetrahedra edges
            fn_CutEdge2Type2(vElem[i], m_vTetrEdge[i], vNewElem);
            break;
        case(Edge2Type3):
            // split two tetrahedra edges
            fn_CutEdge2Type3(vElem[i], m_vTetrEdge[i], vNewElem);
            break;
        case(Edge3Type4):
            // split three tetrahedra edges
            fn_CutEdge3Type4(vElem[i], m_vTetrEdge[i], vNewElem);
            break;
        case(Edge3Type5):
            // split three tetrahedra edges
            fn_CutEdge3Type5(vElem[i], m_vTetrEdge[i], vNewElem);
            break;
        case(Edge3Type6):
            // split three tetrahedra edges
            fn_CutEdge3Type6(vElem[i], m_vTetrEdge[i], vNewElem);
            break;
        case(Edge4Type7):
            // split four tetrahedra edge
            fn_CutEdge4Type7(vElem[i], m_vTetrEdge[i], vNewElem);
            break;
        case(Edge4Type8):
            // split four tetrahedra edge
            fn_CutEdge4Type8(vElem[i], m_vTetrEdge[i], vNewElem);
            break;
        case(Edge5Type9):
            // split five tetrahedra edges
            fn_CutEdge5Type9(vElem[i], m_vTetrEdge[i], vNewElem);
            break;
        case(AllSixEdge):
            // split all six tetrahedra edges
            fn_CutAllSixEdge(vR, vElem[i], m_vTetrEdge[i], vNewElem);
            break;
        default:
            vNewElem.push_back(vElem[i]);
        }
    }
    vNewElem.shrink_to_fit();
    vElem = vNewElem;
}

// split only one tetrahedra edge
void RefineTetr::fn_CutEdge1Type1(Tetr &oldTetr, std::vector<int> &vEdgeIndx, std::vector<Tetr> &vNewTetr)
{
    constexpr int V[6][4] = { {2,0,1,3}, {0,1,2,3}, {1,2,0,3}, {1,0,3,2}, {0,3,1,2}, {0,2,3,1} };
    constexpr int B[6][4] = { {0,3,1,2}, {0,1,2,3}, {0,2,3,1}, {3,2,1,0}, {3,1,0,2}, {2,1,3,0} };
    Tetr tet;

    tet.id = oldTetr.id;

    for (auto i = 0; i < 6; ++i) {
        if (m_vEdgeFlag[vEdgeIndx[i]]) {
            tet.node[0] = oldTetr.node[V[i][0]];
            tet.node[1] = oldTetr.node[V[i][1]];
            tet.node[2] = m_vEdgeNode[vEdgeIndx[i]];
            tet.node[3] = oldTetr.node[V[i][3]];
            tet.boud[0] = oldTetr.boud[B[i][0]];
            tet.boud[1] = oldTetr.boud[B[i][1]];
            tet.boud[2] = -1;
            tet.boud[3] = oldTetr.boud[B[i][3]];
            vNewTetr.push_back(tet);

            tet.node[0] = oldTetr.node[V[i][0]];
            tet.node[1] = m_vEdgeNode[vEdgeIndx[i]];
            tet.node[2] = oldTetr.node[V[i][2]];
            tet.node[3] = oldTetr.node[V[i][3]];
            tet.boud[0] = oldTetr.boud[B[i][0]];
            tet.boud[1] = oldTetr.boud[B[i][1]];
            tet.boud[2] = oldTetr.boud[B[i][2]];
            tet.boud[3] = -1;
            vNewTetr.push_back(tet);
            break;
        }
    }
}

// split two tetrahedra edges
void RefineTetr::fn_CutEdge2Type2(Tetr &oldTetr, std::vector<int> &vEdgeIndx, std::vector<Tetr> &vNewTetr)
{
    int edge[2], node[4], boud[4];
    Tetr tet;

    int m = 0;
    for (auto i = 0; i < 6; ++i) {
        if (m_vEdgeFlag[vEdgeIndx[i]]) {
            edge[m] = i;
            ++m;
        }
    }

    tet.id = oldTetr.id;
    if (edge[0] == 0 && edge[1] == 5) {
        node[0] = oldTetr.node[0];
        node[1] = oldTetr.node[2];
        node[2] = oldTetr.node[3];
        node[3] = oldTetr.node[1];
        boud[0] = oldTetr.boud[2];
        boud[1] = oldTetr.boud[1];
        boud[2] = oldTetr.boud[3];
        boud[3] = oldTetr.boud[0];
        edge[0] = m_vEdgeNode[vEdgeIndx[5]];
        edge[1] = m_vEdgeNode[vEdgeIndx[0]];
    } else if (edge[0] == 1 && edge[1] == 3) {
        node[0] = oldTetr.node[0];
        node[1] = oldTetr.node[1];
        node[2] = oldTetr.node[2];
        node[3] = oldTetr.node[3];
        boud[0] = oldTetr.boud[0];
        boud[1] = oldTetr.boud[1];
        boud[2] = oldTetr.boud[2];
        boud[3] = oldTetr.boud[3];
        edge[0] = m_vEdgeNode[vEdgeIndx[1]];
        edge[1] = m_vEdgeNode[vEdgeIndx[3]];
    } else if (edge[0] == 2 && edge[1] == 4) {
        node[0] = oldTetr.node[0];
        node[1] = oldTetr.node[3];
        node[2] = oldTetr.node[1];
        node[3] = oldTetr.node[2];
        boud[0] = oldTetr.boud[3];
        boud[1] = oldTetr.boud[1];
        boud[2] = oldTetr.boud[0];
        boud[3] = oldTetr.boud[2];
        edge[0] = m_vEdgeNode[vEdgeIndx[4]];
        edge[1] = m_vEdgeNode[vEdgeIndx[2]];
    } else {
        return;
    }
    // deal with type 1
    tet.node[0] = node[0];
    tet.node[1] = node[1];
    tet.node[2] = edge[0];
    tet.node[3] = edge[1];
    tet.boud[0] = boud[0];
    tet.boud[1] = -1;
    tet.boud[2] = -1;
    tet.boud[3] = boud[3];
    vNewTetr.push_back(tet);

    tet.node[0] = node[0];
    tet.node[1] = edge[0];
    tet.node[2] = node[2];
    tet.node[3] = edge[1];
    tet.boud[0] = boud[0];
    tet.boud[1] = -1;
    tet.boud[2] = boud[2];
    tet.boud[3] = -1;
    vNewTetr.push_back(tet);

    tet.node[0] = edge[1];
    tet.node[1] = node[1];
    tet.node[2] = edge[0];
    tet.node[3] = node[3];
    tet.boud[0] = -1;
    tet.boud[1] = boud[1];
    tet.boud[2] = -1;
    tet.boud[3] = boud[3];
    vNewTetr.push_back(tet);

    tet.node[0] = edge[1];
    tet.node[1] = edge[0];
    tet.node[2] = node[2];
    tet.node[3] = node[3];
    tet.boud[0] = -1;
    tet.boud[1] = boud[1];
    tet.boud[2] = boud[2];
    tet.boud[3] = -1;
    vNewTetr.push_back(tet);
}

// split two tetrahedra edges
void RefineTetr::fn_CutEdge2Type3(Tetr &oldTetr, std::vector<int> &vEdgeIndx, std::vector<Tetr> &vNewTetr)
{
    constexpr int EdgeNode[6][2] = { {0,1}, {1,2}, {0,2}, {0,3}, {1,3}, {2,3} };
    constexpr int FaceEdge[4][3] = { {0,1,2}, {1,4,5}, {2,3,5}, {0,3,4} };
    int edge[2], node[4], boud[4];
    bool bType, bEdge1, bEdge2;
    Tetr tet;

    int m = 0;
    for (auto i = 0; i < 6; ++i) {
        if (m_vEdgeFlag[vEdgeIndx[i]]) {
            edge[m] = i;
            ++m;
        }
    }
    boud[0] = EdgeNode[edge[0]][0];
    boud[1] = EdgeNode[edge[0]][1];
    boud[2] = EdgeNode[edge[1]][0];
    boud[3] = EdgeNode[edge[1]][1];
    if (boud[0] == boud[2]) {
        node[0] = oldTetr.node[boud[0]];
        node[1] = oldTetr.node[boud[1]];
        node[2] = oldTetr.node[boud[3]];
    } else if (boud[0] == boud[3]) {
        node[0] = oldTetr.node[boud[0]];
        node[1] = oldTetr.node[boud[1]];
        node[2] = oldTetr.node[boud[2]];
    } else if (boud[1] == boud[2]) {
        node[0] = oldTetr.node[boud[1]];
        node[1] = oldTetr.node[boud[0]];
        node[2] = oldTetr.node[boud[3]];
    } else {
        node[0] = oldTetr.node[boud[1]];
        node[1] = oldTetr.node[boud[0]];
        node[2] = oldTetr.node[boud[2]];
    }
    for (auto i = 0; i < 4; ++i) {
        bType = true;
        for (auto j = 0; j < 3; ++j) {
            if (oldTetr.node[i] == node[j]) {
                bType = false;
                break;
            }
        }
        if (bType) {
            node[3] = oldTetr.node[i];
            break;
        }
    }
    // copy boundary condition
    for (auto i = 0; i < 4; ++i) {
        bEdge1 = false;
        bEdge2 = false;
        for (auto j = 0; j < 3; ++j) {
            if (FaceEdge[i][j] == edge[0]) bEdge1 = true;
            if (FaceEdge[i][j] == edge[1]) bEdge2 = true;
        }
        if (bEdge1 && bEdge2) {
            boud[0] = oldTetr.boud[i];
        } else if (bEdge1) {
            boud[3] = oldTetr.boud[i];
        } else if (bEdge2) {
            boud[2] = oldTetr.boud[i];
        } else {
            boud[1] = oldTetr.boud[i];
        }
    }
    // get edge center node
    edge[0] = vEdgeIndx[edge[0]];
    edge[1] = vEdgeIndx[edge[1]];
    Real d01 = m_vEdgeLen[edge[0]] - m_vEdgeLen[edge[1]];
    if (fabs(d01) < MINGAP) d01 = edge[0] - edge[1];
    edge[0] = m_vEdgeNode[edge[0]];
    edge[1] = m_vEdgeNode[edge[1]];
    // split the mesh
    tet.id = oldTetr.id;
    tet.node[0] = node[0];
    tet.node[1] = edge[0];
    tet.node[2] = edge[1];
    tet.node[3] = node[3];
    tet.boud[0] = boud[0];
    tet.boud[1] = -1;
    tet.boud[2] = boud[2];
    tet.boud[3] = boud[3];
    vNewTetr.push_back(tet);
    if (d01 > 0) {
        tet.node[0] = edge[0];
        tet.node[1] = node[1];
        tet.node[2] = node[2];
        tet.node[3] = node[3];
        tet.boud[0] = boud[0];
        tet.boud[1] = boud[1];
        tet.boud[2] = -1;
        tet.boud[3] = boud[3];
        vNewTetr.push_back(tet);

        tet.node[0] = edge[1];
        tet.node[1] = edge[0];
        tet.node[2] = node[2];
        tet.node[3] = node[3];
        tet.boud[0] = boud[0];
        tet.boud[1] = -1;
        tet.boud[2] = boud[2];
        tet.boud[3] = -1;
        vNewTetr.push_back(tet);
    } else {
        tet.node[0] = edge[0];
        tet.node[1] = node[1];
        tet.node[2] = edge[1];
        tet.node[3] = node[3];
        tet.boud[0] = boud[0];
        tet.boud[1] = -1;
        tet.boud[2] = -1;
        tet.boud[3] = boud[3];
        vNewTetr.push_back(tet);

        tet.node[0] = edge[1];
        tet.node[1] = node[1];
        tet.node[2] = node[2];
        tet.node[3] = node[3];
        tet.boud[0] = boud[0];
        tet.boud[1] = boud[1];
        tet.boud[2] = boud[2];
        tet.boud[3] = -1;
        vNewTetr.push_back(tet);
    }
}

// split three tetrahedra edges
void RefineTetr::fn_CutEdge3Type4(Tetr &oldTetr, std::vector<int> &vEdgeIndx, std::vector<Tetr> &vNewTetr)
{
    int node[4], boud[4], edge[3], e[3], m;
    Tetr tet;

    m = 0;
    for (auto i = 0; i < 6; ++i) {
        if (m_vEdgeFlag[vEdgeIndx[i]]) {
            e[m] = i;
            ++m;
        }
    }
    tet.id = oldTetr.id;
    if (e[0] == 0 && e[1] == 1 && e[2] == 2) {
        node[0] = oldTetr.node[0];
        node[1] = oldTetr.node[1];
        node[2] = oldTetr.node[2];
        node[3] = oldTetr.node[3];
        boud[0] = oldTetr.boud[0];
        boud[1] = oldTetr.boud[1];
        boud[2] = oldTetr.boud[2];
        boud[3] = oldTetr.boud[3];
        edge[0] = m_vEdgeNode[vEdgeIndx[0]];
        edge[1] = m_vEdgeNode[vEdgeIndx[1]];
        edge[2] = m_vEdgeNode[vEdgeIndx[2]];
    } else if (e[0] == 1 && e[1] == 4 && e[2] == 5) {
        node[0] = oldTetr.node[2];
        node[1] = oldTetr.node[1];
        node[2] = oldTetr.node[3];
        node[3] = oldTetr.node[0];
        boud[0] = oldTetr.boud[1];
        boud[1] = oldTetr.boud[3];
        boud[2] = oldTetr.boud[2];
        boud[3] = oldTetr.boud[0];
        edge[0] = m_vEdgeNode[vEdgeIndx[1]];
        edge[1] = m_vEdgeNode[vEdgeIndx[4]];
        edge[2] = m_vEdgeNode[vEdgeIndx[5]];
    } else if (e[0] == 0 && e[1] == 3 && e[2] == 4) {
        node[0] = oldTetr.node[0];
        node[1] = oldTetr.node[3];
        node[2] = oldTetr.node[1];
        node[3] = oldTetr.node[2];
        boud[0] = oldTetr.boud[3];
        boud[1] = oldTetr.boud[1];
        boud[2] = oldTetr.boud[0];
        boud[3] = oldTetr.boud[2];
        edge[0] = m_vEdgeNode[vEdgeIndx[3]];
        edge[1] = m_vEdgeNode[vEdgeIndx[4]];
        edge[2] = m_vEdgeNode[vEdgeIndx[0]];
    } else {
        // e[0] == 2  e[1] == 3  e[2] == 5
        node[0] = oldTetr.node[0];
        node[1] = oldTetr.node[2];
        node[2] = oldTetr.node[3];
        node[3] = oldTetr.node[1];
        boud[0] = oldTetr.boud[2];
        boud[1] = oldTetr.boud[1];
        boud[2] = oldTetr.boud[3];
        boud[3] = oldTetr.boud[0];
        edge[0] = m_vEdgeNode[vEdgeIndx[2]];
        edge[1] = m_vEdgeNode[vEdgeIndx[5]];
        edge[2] = m_vEdgeNode[vEdgeIndx[3]];
    }
    tet.node[0] = node[0];
    tet.node[1] = edge[0];
    tet.node[2] = edge[2];
    tet.node[3] = node[3];
    tet.boud[0] = boud[0];
    tet.boud[1] = -1;
    tet.boud[2] = boud[2];
    tet.boud[3] = boud[3];
    vNewTetr.push_back(tet);

    tet.node[0] = node[1];
    tet.node[1] = edge[1];
    tet.node[2] = edge[0];
    tet.node[3] = node[3];
    tet.boud[0] = boud[0];
    tet.boud[1] = -1;
    tet.boud[2] = boud[3];
    tet.boud[3] = boud[1];
    vNewTetr.push_back(tet);

    tet.node[0] = node[2];
    tet.node[1] = edge[2];
    tet.node[2] = edge[1];
    tet.node[3] = node[3];
    tet.boud[0] = boud[0];
    tet.boud[1] = -1;
    tet.boud[2] = boud[1];
    tet.boud[3] = boud[2];
    vNewTetr.push_back(tet);

    tet.node[0] = edge[0];
    tet.node[1] = edge[1];
    tet.node[2] = edge[2];
    tet.node[3] = node[3];
    tet.boud[0] = boud[0];
    tet.boud[1] = -1;
    tet.boud[2] = -1;
    tet.boud[3] = -1;
    vNewTetr.push_back(tet);
}

// split three tetrahedra edges
void RefineTetr::fn_CutEdge3Type5(Tetr &oldTetr, std::vector<int> &vEdgeIndx, std::vector<Tetr> &vNewTetr)
{
    int edge[3], node[4], boud[4], e[3], m;
    Real d12, delta;
    Tetr tet;

    m = 0;
    for (auto i = 0; i < 6; ++i) {
        if (m_vEdgeFlag[vEdgeIndx[i]]) {
            e[m] = i;
            ++m;
        }
    }
    // find the max edge
    m = e[0];
    delta = m_vEdgeLen[vEdgeIndx[m]];
    for (auto i = 1; i < 3; ++i) {
        d12 = m_vEdgeLen[vEdgeIndx[e[i]]] - delta;
        if (d12 > MINGAP) {
            m = e[i];
            delta = m_vEdgeLen[vEdgeIndx[m]];
        } else if (fabs(d12) < MINGAP) {
            if (vEdgeIndx[e[i]] > vEdgeIndx[m]) {
                m = e[i];
            }
        }
    }
    // get the fourth node
    if (e[0] == 0 && e[1] == 2 && e[2] == 3) {
        node[3] = oldTetr.node[0];
        boud[0] = oldTetr.boud[1];
        edge[0] = vEdgeIndx[m];
        if (m == 0) {           // edge 0 is max
            node[0] = oldTetr.node[1];
            node[1] = oldTetr.node[3];
            node[2] = oldTetr.node[2];
            boud[1] = oldTetr.boud[2];
            boud[2] = oldTetr.boud[0];
            boud[3] = oldTetr.boud[3];
            edge[1] = vEdgeIndx[3];
            edge[2] = vEdgeIndx[2];
        } else if (m == 2) {    // edge 2 is max
            node[0] = oldTetr.node[2];
            node[1] = oldTetr.node[1];
            node[2] = oldTetr.node[3];
            boud[1] = oldTetr.boud[3];
            boud[2] = oldTetr.boud[2];
            boud[3] = oldTetr.boud[0];
            edge[1] = vEdgeIndx[0];
            edge[2] = vEdgeIndx[3];
        } else {                // edge 3 is max
            node[0] = oldTetr.node[3];
            node[1] = oldTetr.node[2];
            node[2] = oldTetr.node[1];
            boud[1] = oldTetr.boud[0];
            boud[2] = oldTetr.boud[3];
            boud[3] = oldTetr.boud[2];
            edge[1] = vEdgeIndx[2];
            edge[2] = vEdgeIndx[0];
        }
    } else if (e[0] == 0 && e[1] == 1 && e[2] == 4) {
        node[3] = oldTetr.node[1];
        boud[0] = oldTetr.boud[2];
        edge[0] = vEdgeIndx[m];
        if (m == 0) {           // edge 0 is max
            node[0] = oldTetr.node[0];
            node[1] = oldTetr.node[2];
            node[2] = oldTetr.node[3];
            boud[1] = oldTetr.boud[1];
            boud[2] = oldTetr.boud[3];
            boud[3] = oldTetr.boud[0];
            edge[1] = vEdgeIndx[1];
            edge[2] = vEdgeIndx[4];
        } else if (m == 1) {    // edge 1 is max
            node[0] = oldTetr.node[2];
            node[1] = oldTetr.node[3];
            node[2] = oldTetr.node[0];
            boud[1] = oldTetr.boud[3];
            boud[2] = oldTetr.boud[0];
            boud[3] = oldTetr.boud[1];
            edge[1] = vEdgeIndx[4];
            edge[2] = vEdgeIndx[0];
        } else {                // edge 4 is max
            node[0] = oldTetr.node[3];
            node[1] = oldTetr.node[0];
            node[2] = oldTetr.node[2];
            boud[1] = oldTetr.boud[0];
            boud[2] = oldTetr.boud[1];
            boud[3] = oldTetr.boud[3];
            edge[1] = vEdgeIndx[0];
            edge[2] = vEdgeIndx[1];
        }
    } else if (e[0] == 1 && e[1] == 2 && e[2] == 5) {
        node[3] = oldTetr.node[2];
        boud[0] = oldTetr.boud[3];
        edge[0] = vEdgeIndx[m];
        if (m == 1) {           // edge 1 is max
            node[0] = oldTetr.node[1];
            node[1] = oldTetr.node[0];
            node[2] = oldTetr.node[3];
            boud[1] = oldTetr.boud[2];
            boud[2] = oldTetr.boud[1];
            boud[3] = oldTetr.boud[0];
            edge[1] = vEdgeIndx[2];
            edge[2] = vEdgeIndx[5];
        } else if (m == 2) {   // edge 2 is max
            node[0] = oldTetr.node[0];
            node[1] = oldTetr.node[3];
            node[2] = oldTetr.node[1];
            boud[1] = oldTetr.boud[1];
            boud[2] = oldTetr.boud[0];
            boud[3] = oldTetr.boud[2];
            edge[1] = vEdgeIndx[5];
            edge[2] = vEdgeIndx[1];
        } else {                // edge 5 is max
            node[0] = oldTetr.node[3];
            node[1] = oldTetr.node[1];
            node[2] = oldTetr.node[0];
            boud[1] = oldTetr.boud[0];
            boud[2] = oldTetr.boud[2];
            boud[3] = oldTetr.boud[1];
            edge[1] = vEdgeIndx[1];
            edge[2] = vEdgeIndx[2];
        }
    } else {
        node[3] = oldTetr.node[3];
        boud[0] = oldTetr.boud[0];
        edge[0] = vEdgeIndx[m];
        if (m == 3) {           // edge 3 is max
            node[0] = oldTetr.node[0];
            node[1] = oldTetr.node[1];
            node[2] = oldTetr.node[2];
            boud[1] = oldTetr.boud[1];
            boud[2] = oldTetr.boud[2];
            boud[3] = oldTetr.boud[3];
            edge[1] = vEdgeIndx[4];
            edge[2] = vEdgeIndx[5];
        } else if (m == 4) {   // edge 4 is max
            node[0] = oldTetr.node[1];
            node[1] = oldTetr.node[2];
            node[2] = oldTetr.node[0];
            boud[1] = oldTetr.boud[2];
            boud[2] = oldTetr.boud[3];
            boud[3] = oldTetr.boud[1];
            edge[1] = vEdgeIndx[5];
            edge[2] = vEdgeIndx[3];
        } else {              // edge 5 is max
            node[0] = oldTetr.node[2];
            node[1] = oldTetr.node[0];
            node[2] = oldTetr.node[1];
            boud[1] = oldTetr.boud[3];
            boud[2] = oldTetr.boud[1];
            boud[3] = oldTetr.boud[2];
            edge[1] = vEdgeIndx[3];
            edge[2] = vEdgeIndx[4];
        }
    }
	delta = m_vEdgeLen[edge[1]] - m_vEdgeLen[edge[2]];
    if (fabs(delta) < MINGAP) delta = edge[1] - edge[2];
	
	edge[0] = m_vEdgeNode[edge[0]];
    edge[1] = m_vEdgeNode[edge[1]];
	edge[2] = m_vEdgeNode[edge[2]];

    tet.id = oldTetr.id;
    tet.node[0] = edge[0]; tet.node[1] = edge[1]; tet.node[2] = edge[2]; tet.node[3] = node[3];
    tet.boud[0] = -1;      tet.boud[1] = boud[1]; tet.boud[2] = boud[2]; tet.boud[3] = boud[3];
    vNewTetr.push_back(tet);

    tet.node[0] = node[0]; tet.node[1] = node[1]; tet.node[2] = node[2]; tet.node[3] = edge[0];
    tet.boud[0] = boud[0]; tet.boud[1] = -1;      tet.boud[2] = boud[2]; tet.boud[3] = boud[3];
    vNewTetr.push_back(tet);

    if (delta > 0) {
        tet.node[0] = node[1]; tet.node[1] = edge[1]; tet.node[2] = node[2]; tet.node[3] = edge[0];
        tet.boud[0] = boud[1]; tet.boud[1] = -1;      tet.boud[2] = -1;      tet.boud[3] = boud[3];
        vNewTetr.push_back(tet);

        tet.node[0] = node[2]; tet.node[1] = edge[1]; tet.node[2] = edge[2]; tet.node[3] = edge[0];
        tet.boud[0] = boud[1]; tet.boud[1] = -1;      tet.boud[2] = boud[2]; tet.boud[3] = -1;
        vNewTetr.push_back(tet);
    } else {
        tet.node[0] = node[1]; tet.node[1] = edge[1]; tet.node[2] = edge[2]; tet.node[3] = edge[0];
        tet.boud[0] = boud[1]; tet.boud[1] = -1;      tet.boud[2] = -1;      tet.boud[3] = boud[3];
        vNewTetr.push_back(tet);

        tet.node[0] = node[1]; tet.node[1] = edge[2]; tet.node[2] = node[2]; tet.node[3] = edge[0];
        tet.boud[0] = boud[1]; tet.boud[1] = boud[2]; tet.boud[2] = -1;      tet.boud[3] = -1;
        vNewTetr.push_back(tet);
    }
}

// split three tetrahedra edges
void RefineTetr::fn_CutEdge3Type6(Tetr &oldTetr, std::vector<int> &vEdgeIndx, std::vector<Tetr> &vNewTetr)
{
    constexpr int EdgeNode[6][2] = { {0,1}, {1,2}, {0,2}, {0,3}, {1,3}, {2,3} };
    constexpr int FaceNode[4][3] = { {0,1,2}, {1,2,3}, {0,2,3}, {0,1,3} };
    int e[3], node[4], boud[4], edge[3], m;
    Real d01, d12;
    Tetr tet;

    boud[0] = 0; boud[1] = 0; boud[2] = 0; boud[3] = 0;
    m = 0;
    for (auto i = 0; i < 6; ++i) {
        if (m_vEdgeFlag[vEdgeIndx[i]]) {
            e[m] = i;
            ++m;
            ++boud[EdgeNode[i][0]];
            ++boud[EdgeNode[i][1]];
        }
    }
    m = 0;
    for (auto i = 0; i < 4; ++i) {
        if (boud[i] == 2) {
            node[m] = i;
            m += 3;
        }
    }
    // get node and edge in new coordinate
    for (auto i = 0; i < 3; ++i) {
        if (EdgeNode[e[i]][0] == node[0] && EdgeNode[e[i]][1] != node[3]) {
            node[1] = EdgeNode[e[i]][1];
            edge[0] = vEdgeIndx[e[i]];
        }
        if (EdgeNode[e[i]][1] == node[0] && EdgeNode[e[i]][0] != node[3]) {
            node[1] = EdgeNode[e[i]][0];
            edge[0] = vEdgeIndx[e[i]];
        }
        if (EdgeNode[e[i]][0] == node[3] && EdgeNode[e[i]][1] != node[0]) {
            node[2] = EdgeNode[e[i]][1];
            edge[2] = vEdgeIndx[e[i]];
        }
        if (EdgeNode[e[i]][1] == node[3] && EdgeNode[e[i]][0] != node[0]) {
            node[2] = EdgeNode[e[i]][0];
            edge[2] = vEdgeIndx[e[i]];
        }
        if (EdgeNode[e[i]][0] == node[0] && EdgeNode[e[i]][1] == node[3]) {
            edge[1] = vEdgeIndx[e[i]];
        }
    }
    d01 = m_vEdgeLen[edge[0]] - m_vEdgeLen[edge[1]];
    d12 = m_vEdgeLen[edge[1]] - m_vEdgeLen[edge[2]];
    if (fabs(d01) < MINGAP) d01 = edge[0] - edge[1];
    if (fabs(d12) < MINGAP) d12 = edge[1] - edge[2];
    // get boundary condition in new coordinate
    for (auto i = 0; i < 4; ++i) {
        m = node[FaceNode[i][0]] + node[FaceNode[i][1]] + node[FaceNode[i][2]];
        if (m == 3) {
            boud[i] = oldTetr.boud[0];
        } else if (m == 6) {
            boud[i] = oldTetr.boud[1];
        } else if (m == 5) {
            boud[i] = oldTetr.boud[2];
        } else {
            boud[i] = oldTetr.boud[3];
        }
    }
    for (auto i = 0; i < 4; ++i) {
        node[i] = oldTetr.node[node[i]];
    }
    edge[0] = m_vEdgeNode[edge[0]];
    edge[1] = m_vEdgeNode[edge[1]];
    edge[2] = m_vEdgeNode[edge[2]];

    // split the element
    tet.id = oldTetr.id;
    tet.node[0] = node[2];
    tet.node[1] = edge[0];
    tet.node[2] = node[1];
    tet.node[3] = edge[2];
    tet.boud[0] = boud[0];
    tet.boud[1] = -1;
    tet.boud[2] = boud[1];
    tet.boud[3] = -1;
    vNewTetr.push_back(tet);

    if (d01 > 0) {
        tet.node[0] = edge[0];
        tet.node[1] = edge[1];
        tet.node[2] = node[3];
        tet.node[3] = edge[2];
        tet.boud[0] = boud[3];
        tet.boud[1] = boud[2];
        tet.boud[2] = -1;
        tet.boud[3] = -1;
        vNewTetr.push_back(tet);

        tet.node[0] = edge[0];
        tet.node[1] = node[3];
        tet.node[2] = node[1];
        tet.node[3] = edge[2];
        tet.boud[0] = boud[3];
        tet.boud[1] = boud[1];
        tet.boud[2] = -1;
        tet.boud[3] = -1;
        vNewTetr.push_back(tet);
    } else {
        tet.node[0] = edge[1];
        tet.node[1] = node[3];
        tet.node[2] = node[1];
        tet.node[3] = edge[2];
        tet.boud[0] = boud[3];
        tet.boud[1] = boud[1];
        tet.boud[2] = -1;
        tet.boud[3] = boud[2];
        vNewTetr.push_back(tet);

        tet.node[0] = edge[1];
        tet.node[1] = node[1];
        tet.node[2] = edge[0];
        tet.node[3] = edge[2];
        tet.boud[0] = boud[3];
        tet.boud[1] = -1;
        tet.boud[2] = -1;
        tet.boud[3] = -1;
        vNewTetr.push_back(tet);
    }

    if (d12 > 0) {
        tet.node[0] = node[0];
        tet.node[1] = node[2];
        tet.node[2] = edge[1];
        tet.node[3] = edge[0];
        tet.boud[0] = boud[2];
        tet.boud[1] = -1;
        tet.boud[2] = boud[3];
        tet.boud[3] = boud[0];
        vNewTetr.push_back(tet);

        tet.node[0] = edge[2];
        tet.node[1] = edge[1];
        tet.node[2] = node[2];
        tet.node[3] = edge[0];
        tet.boud[0] = boud[2];
        tet.boud[1] = -1;
        tet.boud[2] = -1;
        tet.boud[3] = -1;
        vNewTetr.push_back(tet);
    } else {
        tet.node[0] = node[0];
        tet.node[1] = edge[2];
        tet.node[2] = edge[1];
        tet.node[3] = edge[0];
        tet.boud[0] = boud[2];
        tet.boud[1] = -1;
        tet.boud[2] = boud[3];
        tet.boud[3] = -1;
        vNewTetr.push_back(tet);

        tet.node[0] = node[0];
        tet.node[1] = node[2];
        tet.node[2] = edge[2];
        tet.node[3] = edge[0];
        tet.boud[0] = boud[2];
        tet.boud[1] = -1;
        tet.boud[2] = -1;
        tet.boud[3] = boud[0];
        vNewTetr.push_back(tet);
    }
}

// split four tetrahedra edge
void RefineTetr::fn_CutEdge4Type7(Tetr &oldTetr, std::vector<int> &vEdgeIndx, std::vector<Tetr> &vNewTetr)
{
    int e[4], edge[4], node[4], boud[4], m;
    Real delta, d23;
    Tetr tet;

    m = 0;
    for (auto i = 0; i < 6; ++i) {
        if (m_vEdgeFlag[vEdgeIndx[i]]) {
            edge[m] = i;
            ++m;
        }
    }
    m = edge[0];
    delta = m_vEdgeLen[vEdgeIndx[m]];
    for (auto i = 1; i < 4; ++i) {
        d23 = m_vEdgeLen[vEdgeIndx[edge[i]]] - delta;
        if (d23 > MINGAP) {
            m = edge[i];
            delta = m_vEdgeLen[vEdgeIndx[m]];
        } else if (fabs(d23) < MINGAP) {
            if (vEdgeIndx[edge[i]] > vEdgeIndx[m]) {
                m = edge[i];
            }
        }
    }

    if (edge[0] == 1 && edge[1] == 2 && edge[2] == 3 && edge[3] == 4) {
        // split edge 1, 2, 3, 4
        if (m == 1 || m == 2) {
            node[0] = oldTetr.node[2];
            node[1] = oldTetr.node[0];
            node[2] = oldTetr.node[1];
            node[3] = oldTetr.node[3];
            boud[0] = oldTetr.boud[0];
            boud[1] = oldTetr.boud[3];
            boud[2] = oldTetr.boud[1];
            boud[3] = oldTetr.boud[2];
            e[0] = edge[1]; e[1] = edge[0]; e[2] = edge[2]; e[3] = edge[3];
        } else {
            node[0] = oldTetr.node[3];
            node[1] = oldTetr.node[1];
            node[2] = oldTetr.node[0];
            node[3] = oldTetr.node[2];
            boud[0] = oldTetr.boud[3];
            boud[1] = oldTetr.boud[0];
            boud[2] = oldTetr.boud[2];
            boud[3] = oldTetr.boud[1];
            e[0] = edge[3]; e[1] = edge[2]; e[2] = edge[0]; e[3] = edge[1];
        }
    } else if (edge[0] == 0 && edge[1] == 2 && edge[2] == 4 && edge[3] == 5) {
        // split edge 0, 2, 4, 5
        if (m == 0 || m == 2) {
            node[0] = oldTetr.node[0];
            node[1] = oldTetr.node[1];
            node[2] = oldTetr.node[2];
            node[3] = oldTetr.node[3];
            boud[0] = oldTetr.boud[0];
            boud[1] = oldTetr.boud[1];
            boud[2] = oldTetr.boud[2];
            boud[3] = oldTetr.boud[3];
            e[0] = edge[0]; e[1] = edge[1]; e[2] = edge[2]; e[3] = edge[3];
        } else {
            node[0] = oldTetr.node[3];
            node[1] = oldTetr.node[2];
            node[2] = oldTetr.node[1];
            node[3] = oldTetr.node[0];
            boud[0] = oldTetr.boud[1];
            boud[1] = oldTetr.boud[0];
            boud[2] = oldTetr.boud[3];
            boud[3] = oldTetr.boud[2];
            e[0] = edge[3]; e[1] = edge[2]; e[2] = edge[1]; e[3] = edge[0];
        }
    } else {
        // split edge 0, 1, 3, 5
        if (m == 0 || m == 1) {
            node[0] = oldTetr.node[1];
            node[1] = oldTetr.node[2];
            node[2] = oldTetr.node[0];
            node[3] = oldTetr.node[3];
            boud[0] = oldTetr.boud[0]; 
            boud[1] = oldTetr.boud[2];
            boud[2] = oldTetr.boud[3];
            boud[3] = oldTetr.boud[1];
            e[0] = edge[1]; e[1] = edge[0]; e[2] = edge[3]; e[3] = edge[2];
        } else {
            node[0] = oldTetr.node[3];
            node[1] = oldTetr.node[0];
            node[2] = oldTetr.node[2];
            node[3] = oldTetr.node[1];
            boud[0] = oldTetr.boud[2];
            boud[1] = oldTetr.boud[0];
            boud[2] = oldTetr.boud[1];
            boud[3] = oldTetr.boud[3];
            e[0] = edge[2]; e[1] = edge[3]; e[2] = edge[0]; e[3] = edge[1];
        }
    }
    edge[0] = m_vEdgeNode[vEdgeIndx[e[0]]];
    edge[1] = m_vEdgeNode[vEdgeIndx[e[1]]];
    edge[2] = m_vEdgeNode[vEdgeIndx[e[2]]];
    edge[3] = m_vEdgeNode[vEdgeIndx[e[3]]];
    // split the mesh
    tet.id = oldTetr.id;

    d23 = m_vEdgeLen[vEdgeIndx[e[2]]] - m_vEdgeLen[vEdgeIndx[e[3]]];
    if (fabs(d23) < MINGAP) d23 = vEdgeIndx[e[2]] - vEdgeIndx[e[3]];

    if (m == e[0]) {
        delta = m_vEdgeLen[vEdgeIndx[e[1]]] - m_vEdgeLen[vEdgeIndx[e[3]]];
        if (fabs(delta) < MINGAP) delta = vEdgeIndx[e[1]] - vEdgeIndx[e[3]];

        tet.node[0] = edge[1]; tet.node[1] = edge[0]; tet.node[2] = node[2]; tet.node[3] = edge[3];
        tet.boud[0] = boud[0]; tet.boud[1] = -1; tet.boud[2] = boud[2]; tet.boud[3] = -1;
        vNewTetr.push_back(tet);

        tet.node[0] = edge[0]; tet.node[1] = edge[3]; tet.node[2] = node[3]; tet.node[3] = edge[2];
        tet.boud[0] = -1; tet.boud[1] = boud[1]; tet.boud[2] = boud[3]; tet.boud[3] = -1;
        vNewTetr.push_back(tet);

        if (delta > 0) {
            tet.node[0] = node[0]; tet.node[1] = edge[0]; tet.node[2] = edge[1]; tet.node[3] = node[3];
            tet.boud[0] = boud[0]; tet.boud[1] = -1; tet.boud[2] = boud[2]; tet.boud[3] = boud[3];
            vNewTetr.push_back(tet);

            tet.node[0] = edge[3]; tet.node[1] = edge[1]; tet.node[2] = edge[0]; tet.node[3] = node[3];
            tet.boud[0] = -1; tet.boud[1] = -1; tet.boud[2] = -1; tet.boud[3] = boud[2];
            vNewTetr.push_back(tet);
        } else {
            tet.node[0] = node[0]; tet.node[1] = edge[1]; tet.node[2] = edge[3]; tet.node[3] = edge[0];
            tet.boud[0] = boud[2]; tet.boud[1] = -1; tet.boud[2] = -1; tet.boud[3] = boud[0];
            vNewTetr.push_back(tet);

            tet.node[0] = node[0]; tet.node[1] = edge[3]; tet.node[2] = node[3]; tet.node[3] = edge[0];
            tet.boud[0] = boud[2]; tet.boud[1] = -1; tet.boud[2] = boud[3]; tet.boud[3] = -1;
            vNewTetr.push_back(tet);
        }

        if (d23 > 0) {
            tet.node[0] = edge[0]; tet.node[1] = node[1]; tet.node[2] = node[2]; tet.node[3] = edge[2];
            tet.boud[0] = boud[0]; tet.boud[1] = boud[1]; tet.boud[2] = -1; tet.boud[3] = boud[3];
            vNewTetr.push_back(tet);

            tet.node[0] = edge[0]; tet.node[1] = node[2]; tet.node[2] = edge[3]; tet.node[3] = edge[2];
            tet.boud[0] = -1; tet.boud[1] = boud[1]; tet.boud[2] = -1; tet.boud[3] = -1;
            vNewTetr.push_back(tet);
        } else {
            tet.node[0] = edge[0]; tet.node[1] = node[1]; tet.node[2] = node[2]; tet.node[3] = edge[3];
            tet.boud[0] = boud[0]; tet.boud[1] = boud[1]; tet.boud[2] = -1; tet.boud[3] = -1;
            vNewTetr.push_back(tet);

            tet.node[0] = edge[0]; tet.node[1] = node[1]; tet.node[2] = edge[3]; tet.node[3] = edge[2];
            tet.boud[0] = -1; tet.boud[1] = boud[1]; tet.boud[2] = -1; tet.boud[3] = boud[3];
            vNewTetr.push_back(tet);
        }
    } else {
        delta = m_vEdgeLen[vEdgeIndx[e[0]]] - m_vEdgeLen[vEdgeIndx[e[2]]];
        if (fabs(delta) < MINGAP) delta = vEdgeIndx[e[0]] - vEdgeIndx[e[2]];

        tet.node[0] = edge[0]; tet.node[1] = node[1]; tet.node[2] = edge[1]; tet.node[3] = edge[2];
        tet.boud[0] = boud[0]; tet.boud[1] = -1; tet.boud[2] = -1; tet.boud[3] = boud[3];
        vNewTetr.push_back(tet);

        tet.node[0] = edge[1]; tet.node[1] = edge[2]; tet.node[2] = edge[3]; tet.node[3] = node[3];
        tet.boud[0] = -1; tet.boud[1] = boud[1]; tet.boud[2] = boud[2]; tet.boud[3] = -1;
        vNewTetr.push_back(tet);

        if (delta > 0) {
            tet.node[0] = edge[2]; tet.node[1] = edge[0]; tet.node[2] = node[3]; tet.node[3] = edge[1];
            tet.boud[0] = boud[3]; tet.boud[1] = -1; tet.boud[2] = -1; tet.boud[3] = -1;
            vNewTetr.push_back(tet);

            tet.node[0] = node[0]; tet.node[1] = node[3]; tet.node[2] = edge[0]; tet.node[3] = edge[1];
            tet.boud[0] = boud[3]; tet.boud[1] = -1; tet.boud[2] = boud[0]; tet.boud[3] = boud[2];
            vNewTetr.push_back(tet);
        } else {
            tet.node[0] = node[0]; tet.node[1] = edge[2]; tet.node[2] = edge[0]; tet.node[3] = edge[1];
            tet.boud[0] = boud[3]; tet.boud[1] = -1; tet.boud[2] = boud[0]; tet.boud[3] = -1;
            vNewTetr.push_back(tet);

            tet.node[0] = node[0]; tet.node[1] = node[3]; tet.node[2] = edge[2]; tet.node[3] = edge[1];
            tet.boud[0] = boud[3]; tet.boud[1] = -1; tet.boud[2] = -1; tet.boud[3] = boud[2];
            vNewTetr.push_back(tet);
        }

        if (d23 > 0) {
            tet.node[0] = edge[1]; tet.node[1] = edge[3]; tet.node[2] = edge[2]; tet.node[3] = node[2];
            tet.boud[0] = -1; tet.boud[1] = boud[1]; tet.boud[2] = -1; tet.boud[3] = boud[2];
            vNewTetr.push_back(tet);

            tet.node[0] = edge[1]; tet.node[1] = edge[2]; tet.node[2] = node[1]; tet.node[3] = node[2];
            tet.boud[0] = -1; tet.boud[1] = boud[1]; tet.boud[2] = boud[0]; tet.boud[3] = -1;
            vNewTetr.push_back(tet);
        } else {
            tet.node[0] = edge[1]; tet.node[1] = node[2]; tet.node[2] = edge[3]; tet.node[3] = node[1];
            tet.boud[0] = boud[2]; tet.boud[1] = boud[1]; tet.boud[2] = -1; tet.boud[3] = boud[0];
            vNewTetr.push_back(tet);

            tet.node[0] = edge[1]; tet.node[1] = edge[3]; tet.node[2] = edge[2]; tet.node[3] = node[1];
            tet.boud[0] = -1; tet.boud[1] = boud[1]; tet.boud[2] = -1; tet.boud[3] = -1;
            vNewTetr.push_back(tet);
        }
    }
}

// split four tetrahedra edge
void RefineTetr::fn_CutEdge4Type8(Tetr &oldTetr, std::vector<int> &vEdgeIndx, std::vector<Tetr> &vNewTetr)
{
    int edge[4], node[4], boud[4];
    Real d03, d23;
    Tetr tet;

    if (m_vEdgeFlag[vEdgeIndx[0]] && m_vEdgeFlag[vEdgeIndx[1]] && m_vEdgeFlag[vEdgeIndx[2]]) {
        node[3] = oldTetr.node[3];
        boud[0] = oldTetr.boud[0];
        if (m_vEdgeFlag[vEdgeIndx[3]]) {
            node[0] = oldTetr.node[0]; node[1] = oldTetr.node[1]; node[2] = oldTetr.node[2];
            boud[1] = oldTetr.boud[1]; boud[2] = oldTetr.boud[2]; boud[3] = oldTetr.boud[3];
            edge[0] = vEdgeIndx[0];
            edge[1] = vEdgeIndx[1];
            edge[2] = vEdgeIndx[2];
            edge[3] = vEdgeIndx[3];
        } else if (m_vEdgeFlag[vEdgeIndx[4]]) {
            node[0] = oldTetr.node[1]; node[1] = oldTetr.node[2]; node[2] = oldTetr.node[0];
            boud[1] = oldTetr.boud[2]; boud[2] = oldTetr.boud[3]; boud[3] = oldTetr.boud[1];
            edge[0] = vEdgeIndx[1];
            edge[1] = vEdgeIndx[2];
            edge[2] = vEdgeIndx[0];
            edge[3] = vEdgeIndx[4];
        } else {
            node[0] = oldTetr.node[2]; node[1] = oldTetr.node[0]; node[2] = oldTetr.node[1];
            boud[1] = oldTetr.boud[3]; boud[2] = oldTetr.boud[1]; boud[3] = oldTetr.boud[2];
            edge[0] = vEdgeIndx[2];
            edge[1] = vEdgeIndx[0];
            edge[2] = vEdgeIndx[1];
            edge[3] = vEdgeIndx[5];
        }
    } else if (m_vEdgeFlag[vEdgeIndx[1]] && m_vEdgeFlag[vEdgeIndx[4]] && m_vEdgeFlag[vEdgeIndx[5]]) {
        node[3] = oldTetr.node[0];
        boud[0] = oldTetr.boud[1];
        if (m_vEdgeFlag[vEdgeIndx[0]]) {
            node[0] = oldTetr.node[1]; node[1] = oldTetr.node[3]; node[2] = oldTetr.node[2];
            boud[1] = oldTetr.boud[2]; boud[2] = oldTetr.boud[0]; boud[3] = oldTetr.boud[3];
            edge[0] = vEdgeIndx[4];
            edge[1] = vEdgeIndx[5];
            edge[2] = vEdgeIndx[1];
            edge[3] = vEdgeIndx[0];
        } else if (m_vEdgeFlag[vEdgeIndx[2]]) {
            node[0] = oldTetr.node[2]; node[1] = oldTetr.node[1]; node[2] = oldTetr.node[3];
            boud[1] = oldTetr.boud[3]; boud[2] = oldTetr.boud[2]; boud[3] = oldTetr.boud[0];
            edge[0] = vEdgeIndx[1];
            edge[1] = vEdgeIndx[4]; 
            edge[2] = vEdgeIndx[5];
            edge[3] = vEdgeIndx[2];
        } else {
            node[0] = oldTetr.node[3]; node[1] = oldTetr.node[2]; node[2] = oldTetr.node[1];
            boud[1] = oldTetr.boud[0]; boud[2] = oldTetr.boud[3]; boud[3] = oldTetr.boud[2];
            edge[0] = vEdgeIndx[5];
            edge[1] = vEdgeIndx[1];
            edge[2] = vEdgeIndx[4];
            edge[3] = vEdgeIndx[3];
        }
    } else if(m_vEdgeFlag[vEdgeIndx[2]] && m_vEdgeFlag[vEdgeIndx[3]] && m_vEdgeFlag[vEdgeIndx[5]]) {
        node[3] = oldTetr.node[1];
        boud[0] = oldTetr.boud[2];
        if (m_vEdgeFlag[vEdgeIndx[0]]) {
            node[0] = oldTetr.node[0]; node[1] = oldTetr.node[2]; node[2] = oldTetr.node[3];
            boud[1] = oldTetr.boud[1]; boud[2] = oldTetr.boud[3]; boud[3] = oldTetr.boud[0];
            edge[0] = vEdgeIndx[2];
            edge[1] = vEdgeIndx[5];
            edge[2] = vEdgeIndx[3];
            edge[3] = vEdgeIndx[0];
        } else if (m_vEdgeFlag[vEdgeIndx[1]]) {
            node[0] = oldTetr.node[2]; node[1] = oldTetr.node[3]; node[2] = oldTetr.node[0];
            boud[1] = oldTetr.boud[3]; boud[2] = oldTetr.boud[0]; boud[3] = oldTetr.boud[1];
            edge[0] = vEdgeIndx[5];
            edge[1] = vEdgeIndx[3];
            edge[2] = vEdgeIndx[2];
            edge[3] = vEdgeIndx[1];
        } else {
            node[0] = oldTetr.node[3]; node[1] = oldTetr.node[0]; node[2] = oldTetr.node[2];
            boud[1] = oldTetr.boud[0]; boud[2] = oldTetr.boud[1]; boud[3] = oldTetr.boud[3];
            edge[0] = vEdgeIndx[3];
            edge[1] = vEdgeIndx[2];
            edge[2] = vEdgeIndx[5];
            edge[3] = vEdgeIndx[4];
        }
    } else {
        node[3] = oldTetr.node[2];
        boud[0] = oldTetr.boud[3];
        if (m_vEdgeFlag[vEdgeIndx[1]]) {
            node[0] = oldTetr.node[1]; node[1] = oldTetr.node[0]; node[2] = oldTetr.node[3];
            boud[1] = oldTetr.boud[2]; boud[2] = oldTetr.boud[1]; boud[3] = oldTetr.boud[0];
            edge[0] = vEdgeIndx[0];
            edge[1] = vEdgeIndx[3];
            edge[2] = vEdgeIndx[4]; 
            edge[3] = vEdgeIndx[1];
        } else if (m_vEdgeFlag[vEdgeIndx[2]]) {
            node[0] = oldTetr.node[0]; node[1] = oldTetr.node[3]; node[2] = oldTetr.node[1];
            boud[1] = oldTetr.boud[1]; boud[2] = oldTetr.boud[0]; boud[3] = oldTetr.boud[2];
            edge[0] = vEdgeIndx[3];
            edge[1] = vEdgeIndx[4];
            edge[2] = vEdgeIndx[0];
            edge[3] = vEdgeIndx[2];
        } else {
            node[0] = oldTetr.node[3]; node[1] = oldTetr.node[1]; node[2] = oldTetr.node[0];
            boud[1] = oldTetr.boud[0]; boud[2] = oldTetr.boud[2]; boud[3] = oldTetr.boud[1];
            edge[0] = vEdgeIndx[4];
            edge[1] = vEdgeIndx[0];
            edge[2] = vEdgeIndx[3];
            edge[3] = vEdgeIndx[5];
        }
    }

    d03 = m_vEdgeLen[edge[0]] - m_vEdgeLen[edge[3]];
    d23 = m_vEdgeLen[edge[2]] - m_vEdgeLen[edge[3]];
    if (fabs(d03) < MINGAP) d03 = edge[0] - edge[3];
    if (fabs(d23) < MINGAP) d23 = edge[2] - edge[3];
    edge[0] = m_vEdgeNode[edge[0]];
    edge[1] = m_vEdgeNode[edge[1]];
    edge[2] = m_vEdgeNode[edge[2]];
    edge[3] = m_vEdgeNode[edge[3]];

    // split the element
    tet.id = oldTetr.id;
    tet.node[0] = node[0];
    tet.node[1] = edge[0];
    tet.node[2] = edge[2];
    tet.node[3] = edge[3];
    tet.boud[0] = boud[0];
    tet.boud[1] = -1;
    tet.boud[2] = boud[2];
    tet.boud[3] = boud[3];
    vNewTetr.push_back(tet);

    tet.node[0] = edge[1];
    tet.node[1] = edge[2];
    tet.node[2] = edge[0];
    tet.node[3] = edge[3];
    tet.boud[0] = boud[0];
    tet.boud[1] = -1;
    tet.boud[2] = -1;
    tet.boud[3] = -1;
    vNewTetr.push_back(tet);

    if (d03 > 0) {
        tet.node[0] = edge[3];
        tet.node[1] = edge[0];
        tet.node[2] = edge[1];
        tet.node[3] = node[3];
        tet.boud[0] = -1;
        tet.boud[1] = -1;
        tet.boud[2] = -1;
        tet.boud[3] = boud[3];
        vNewTetr.push_back(tet);

        tet.node[0] = edge[0];
        tet.node[1] = node[1];
        tet.node[2] = edge[1];
        tet.node[3] = node[3];
        tet.boud[0] = boud[0];
        tet.boud[1] = boud[1];
        tet.boud[2] = -1;
        tet.boud[3] = boud[3];
        vNewTetr.push_back(tet);
    } else {
        tet.node[0] = edge[3];
        tet.node[1] = node[1];
        tet.node[2] = edge[1];
        tet.node[3] = node[3];
        tet.boud[0] = -1;
        tet.boud[1] = boud[1];
        tet.boud[2] = -1;
        tet.boud[3] = boud[3];
        vNewTetr.push_back(tet);

        tet.node[0] = edge[0];
        tet.node[1] = edge[3];
        tet.node[2] = node[1];
        tet.node[3] = edge[1];
        tet.boud[0] = boud[3];
        tet.boud[1] = -1;
        tet.boud[2] = boud[0];
        tet.boud[3] = -1;
        vNewTetr.push_back(tet);
    }

    if (d23 > 0) {
        tet.node[0] = edge[3];
        tet.node[1] = edge[1];
        tet.node[2] = edge[2];
        tet.node[3] = node[3];
        tet.boud[0] = -1;
        tet.boud[1] = -1;
        tet.boud[2] = boud[2];
        tet.boud[3] = -1;
        vNewTetr.push_back(tet);

        tet.node[0] = edge[2];
        tet.node[1] = node[2];
        tet.node[2] = node[3];
        tet.node[3] = edge[1];
        tet.boud[0] = boud[2];
        tet.boud[1] = boud[1];
        tet.boud[2] = -1;
        tet.boud[3] = boud[0];
        vNewTetr.push_back(tet);
    } else {
        tet.node[0] = edge[3];
        tet.node[1] = edge[1];
        tet.node[2] = node[2];
        tet.node[3] = node[3];
        tet.boud[0] = -1;
        tet.boud[1] = boud[1];
        tet.boud[2] = boud[2];
        tet.boud[3] = -1;
        vNewTetr.push_back(tet);

        tet.node[0] = edge[2];
        tet.node[1] = node[2];
        tet.node[2] = edge[3];
        tet.node[3] = edge[1];
        tet.boud[0] = boud[2];
        tet.boud[1] = -1;
        tet.boud[2] = -1;
        tet.boud[3] = boud[0];
        vNewTetr.push_back(tet);
    }
}

// split five tetrahedra edges
void RefineTetr::fn_CutEdge5Type9(Tetr &oldTetr, std::vector<int> &vEdgeIndx, std::vector<Tetr> &vNewTetr)
{
    constexpr int V[6][4] = { {2,0,1,3}, {0,1,2,3}, {1,2,0,3}, {1,0,3,2}, {0,3,1,2}, {0,2,3,1} };
    constexpr int B[6][4] = { {0,3,1,2}, {0,1,2,3}, {0,2,3,1}, {3,2,1,0}, {3,1,0,2}, {2,1,3,0} };
    constexpr int E[6][6] = { {2,0,1,5,3,4}, {0,1,2,3,4,5}, {1,2,0,4,5,3}, {0,3,4,1,2,5}, {3,4,0,2,5,1}, {2,5,3,0,1,4} };
    int node[4], boud[4], edge[6];
    Real d02 = 0, d45 = 0;
    Tetr tet;

    for (auto i = 0; i < 6; ++i) {
        if (!m_vEdgeFlag[vEdgeIndx[i]]) {
            node[0] = oldTetr.node[V[i][0]];
            node[1] = oldTetr.node[V[i][1]];
            node[2] = oldTetr.node[V[i][2]];
            node[3] = oldTetr.node[V[i][3]];
            boud[0] = oldTetr.boud[B[i][0]];
            boud[1] = oldTetr.boud[B[i][1]];
            boud[2] = oldTetr.boud[B[i][2]];
            boud[3] = oldTetr.boud[B[i][3]];
            edge[0] = m_vEdgeNode[vEdgeIndx[E[i][0]]];
            edge[1] = m_vEdgeNode[vEdgeIndx[E[i][1]]];
            edge[2] = m_vEdgeNode[vEdgeIndx[E[i][2]]];
            edge[3] = m_vEdgeNode[vEdgeIndx[E[i][3]]];
            edge[4] = m_vEdgeNode[vEdgeIndx[E[i][4]]];
            edge[5] = m_vEdgeNode[vEdgeIndx[E[i][5]]];
            d02 = m_vEdgeLen[vEdgeIndx[E[i][0]]] - m_vEdgeLen[vEdgeIndx[E[i][2]]];
            d45 = m_vEdgeLen[vEdgeIndx[E[i][4]]] - m_vEdgeLen[vEdgeIndx[E[i][5]]];
            if (fabs(d02) < MINGAP) d02 = vEdgeIndx[E[i][0]] - vEdgeIndx[E[i][2]];
            if (fabs(d45) < MINGAP) d45 = vEdgeIndx[E[i][4]] - vEdgeIndx[E[i][5]];
            break;
        }
    }

    tet.id = oldTetr.id;
    tet.node[0] = edge[3];
    tet.node[1] = edge[4];
    tet.node[2] = edge[5];
    tet.node[3] = node[3];
    tet.boud[0] = -1;
    tet.boud[1] = boud[1];
    tet.boud[2] = boud[2];
    tet.boud[3] = boud[3];
    vNewTetr.push_back(tet);

    tet.node[0] = node[0];
    tet.node[1] = edge[0];
    tet.node[2] = edge[2];
    tet.node[3] = edge[3];
    tet.boud[0] = boud[0];
    tet.boud[1] = -1;
    tet.boud[2] = boud[2];
    tet.boud[3] = boud[3];
    vNewTetr.push_back(tet);

    if (d02 > 0) {
        tet.node[0] = edge[0];
        tet.node[1] = edge[5];
        tet.node[2] = edge[3];
        tet.node[3] = edge[4];
        tet.boud[0] = -1;
        tet.boud[1] = -1;
        tet.boud[2] = boud[3];
        tet.boud[3] = -1;
        vNewTetr.push_back(tet);

        tet.node[0] = edge[0];
        tet.node[1] = edge[2];
        tet.node[2] = edge[3];
        tet.node[3] = edge[5];
        tet.boud[0] = -1;
        tet.boud[1] = boud[2];
        tet.boud[2] = -1;
        tet.boud[3] = -1;
        vNewTetr.push_back(tet);

        tet.node[0] = edge[2];
        tet.node[1] = edge[0];
        tet.node[2] = node[2];
        tet.node[3] = edge[5];
        tet.boud[0] = boud[0];
        tet.boud[1] = -1;
        tet.boud[2] = boud[2];
        tet.boud[3] = -1;
        vNewTetr.push_back(tet);
        if (d45 > 0) {
            tet.node[0] = node[1]; tet.node[1] = node[2]; tet.node[2] = edge[0]; tet.node[3] = edge[4];
            tet.boud[0] = boud[0]; tet.boud[1] = -1;      tet.boud[2] = boud[3]; tet.boud[3] = boud[1];
            vNewTetr.push_back(tet);

            tet.node[0] = edge[0]; tet.node[1] = edge[5]; tet.node[2] = edge[4]; tet.node[3] = node[2];
            tet.boud[0] = -1;      tet.boud[1] = boud[1]; tet.boud[2] = -1;      tet.boud[3] = -1;
            vNewTetr.push_back(tet);
        } else {
            tet.node[0] = edge[0]; tet.node[1] = edge[5]; tet.node[2] = edge[4]; tet.node[3] = node[1];
            tet.boud[0] = -1;      tet.boud[1] = boud[1]; tet.boud[2] = boud[3];      tet.boud[3] = -1;
            vNewTetr.push_back(tet);

            tet.node[0] = edge[0]; tet.node[1] = node[2]; tet.node[2] = edge[5]; tet.node[3] = node[1];
            tet.boud[0] = -1;      tet.boud[1] = boud[1]; tet.boud[2] = -1;      tet.boud[3] = boud[0];
            vNewTetr.push_back(tet);
        }
    } else {
        tet.node[0] = edge[0];
        tet.node[1] = edge[2];
        tet.node[2] = edge[3];
        tet.node[3] = edge[4];
        tet.boud[0] = -1;
        tet.boud[1] = -1;
        tet.boud[2] = boud[3];
        tet.boud[3] = -1;
        vNewTetr.push_back(tet);

        tet.node[0] = edge[3];
        tet.node[1] = edge[4];
        tet.node[2] = edge[2];
        tet.node[3] = edge[5];
        tet.boud[0] = -1;
        tet.boud[1] = -1;
        tet.boud[2] = boud[2];
        tet.boud[3] = -1;
        vNewTetr.push_back(tet);

        tet.node[0] = edge[0];
        tet.node[1] = node[1];
        tet.node[2] = edge[2];
        tet.node[3] = edge[4];
        tet.boud[0] = boud[0];
        tet.boud[1] = -1;
        tet.boud[2] = -1;
        tet.boud[3] = boud[3];
        vNewTetr.push_back(tet);
        if (d45 > 0) {
            tet.node[0] = edge[2]; tet.node[1] = edge[4]; tet.node[2] = node[1]; tet.node[3] = node[2];
            tet.boud[0] = -1;      tet.boud[1] = boud[1]; tet.boud[2] = boud[0];      tet.boud[3] = -1;
            vNewTetr.push_back(tet);

            tet.node[0] = edge[2]; tet.node[1] = edge[5]; tet.node[2] = edge[4]; tet.node[3] = node[2];
            tet.boud[0] = -1;      tet.boud[1] = boud[1]; tet.boud[2] = -1;      tet.boud[3] = boud[2];
            vNewTetr.push_back(tet);
        } else {
            tet.node[0] = edge[5]; tet.node[1] = edge[4]; tet.node[2] = edge[2]; tet.node[3] = node[1];
            tet.boud[0] = -1;      tet.boud[1] = -1;      tet.boud[2] = -1;      tet.boud[3] = boud[1];
            vNewTetr.push_back(tet);

            tet.node[0] = edge[2]; tet.node[1] = node[1]; tet.node[2] = node[2]; tet.node[3] = edge[5];
            tet.boud[0] = boud[0]; tet.boud[1] = boud[1]; tet.boud[2] = boud[2];      tet.boud[3] = -1;
            vNewTetr.push_back(tet);
        }
    }
}

// split all six tetrahedra edges
void RefineTetr::fn_CutAllSixEdge(std::vector<Node> &vR, Tetr &oldTetr, std::vector<int> &vEdgeIndx, std::vector<Tetr> &vNewTetr)
{
    Real d05, d13;
    Tetr tet;

    tet.id = oldTetr.id;
    tet.node[0] = m_vEdgeNode[vEdgeIndx[3]];
    tet.node[1] = m_vEdgeNode[vEdgeIndx[4]];
    tet.node[2] = m_vEdgeNode[vEdgeIndx[5]];
    tet.node[3] = oldTetr.node[3];
    tet.boud[0] = -1;
    tet.boud[1] = oldTetr.boud[1];
    tet.boud[2] = oldTetr.boud[2];
    tet.boud[3] = oldTetr.boud[3];
    vNewTetr.push_back(tet);

    tet.node[0] = m_vEdgeNode[vEdgeIndx[0]];
    tet.node[1] = oldTetr.node[1];
    tet.node[2] = m_vEdgeNode[vEdgeIndx[1]];
    tet.node[3] = m_vEdgeNode[vEdgeIndx[4]];
    tet.boud[0] = oldTetr.boud[0];
    tet.boud[1] = oldTetr.boud[1];
    tet.boud[2] = -1;
    tet.boud[3] = oldTetr.boud[3];
    vNewTetr.push_back(tet);

    tet.node[0] = m_vEdgeNode[vEdgeIndx[2]];
    tet.node[1] = m_vEdgeNode[vEdgeIndx[1]];
    tet.node[2] = oldTetr.node[2];
    tet.node[3] = m_vEdgeNode[vEdgeIndx[5]];
    tet.boud[0] = oldTetr.boud[0];
    tet.boud[1] = oldTetr.boud[1];
    tet.boud[2] = oldTetr.boud[2];
    tet.boud[3] = -1;
    vNewTetr.push_back(tet);

    tet.node[0] = oldTetr.node[0];
    tet.node[1] = m_vEdgeNode[vEdgeIndx[0]];
    tet.node[2] = m_vEdgeNode[vEdgeIndx[2]];
    tet.node[3] = m_vEdgeNode[vEdgeIndx[3]];
    tet.boud[0] = oldTetr.boud[0];
    tet.boud[1] = -1;
    tet.boud[2] = oldTetr.boud[2];
    tet.boud[3] = oldTetr.boud[3];
    vNewTetr.push_back(tet);

    Node rd;
    int m = m_vEdgeNode[vEdgeIndx[0]];
    int n = m_vEdgeNode[vEdgeIndx[5]];
    rd.x = vR[m].x - vR[n].x;
    rd.y = vR[m].y - vR[n].y;
    rd.z = vR[m].z - vR[n].z;
    d05 = rd.x * rd.x + rd.y * rd.y + rd.z * rd.z;

    m = m_vEdgeNode[vEdgeIndx[1]];
    n = m_vEdgeNode[vEdgeIndx[3]];
    rd.x = vR[m].x - vR[n].x;
    rd.y = vR[m].y - vR[n].y;
    rd.z = vR[m].z - vR[n].z;
    d13 = rd.x * rd.x + rd.y * rd.y + rd.z * rd.z;

    if (d05 > d13) {
        tet.node[0] = m_vEdgeNode[vEdgeIndx[0]];
        tet.node[1] = m_vEdgeNode[vEdgeIndx[1]];
        tet.node[2] = m_vEdgeNode[vEdgeIndx[3]];
        tet.node[3] = m_vEdgeNode[vEdgeIndx[4]];
        tet.boud[0] = -1;
        tet.boud[1] = -1;
        tet.boud[2] = oldTetr.boud[3];
        tet.boud[3] = -1;
        vNewTetr.push_back(tet);

        tet.node[0] = m_vEdgeNode[vEdgeIndx[5]];
        tet.node[1] = m_vEdgeNode[vEdgeIndx[3]];
        tet.node[2] = m_vEdgeNode[vEdgeIndx[1]];
        tet.node[3] = m_vEdgeNode[vEdgeIndx[4]];
        tet.boud[0] = -1;
        tet.boud[1] = -1;
        tet.boud[2] = oldTetr.boud[1];
        tet.boud[3] = -1;
        vNewTetr.push_back(tet);

        tet.node[0] = m_vEdgeNode[vEdgeIndx[1]];
        tet.node[1] = m_vEdgeNode[vEdgeIndx[2]];
        tet.node[2] = m_vEdgeNode[vEdgeIndx[0]];
        tet.node[3] = m_vEdgeNode[vEdgeIndx[3]];
        tet.boud[0] = oldTetr.boud[0];
        tet.boud[1] = -1;
        tet.boud[2] = -1;
        tet.boud[3] = -1;
        vNewTetr.push_back(tet);

        tet.node[0] = m_vEdgeNode[vEdgeIndx[3]];
        tet.node[1] = m_vEdgeNode[vEdgeIndx[1]];
        tet.node[2] = m_vEdgeNode[vEdgeIndx[2]];
        tet.node[3] = m_vEdgeNode[vEdgeIndx[5]];
        tet.boud[0] = -1;
        tet.boud[1] = -1;
        tet.boud[2] = oldTetr.boud[2];
        tet.boud[3] = -1;
        vNewTetr.push_back(tet);
    } else {
        tet.node[0] = m_vEdgeNode[vEdgeIndx[3]];
        tet.node[1] = m_vEdgeNode[vEdgeIndx[0]];
        tet.node[2] = m_vEdgeNode[vEdgeIndx[5]];
        tet.node[3] = m_vEdgeNode[vEdgeIndx[4]];
        tet.boud[0] = -1;
        tet.boud[1] = -1;
        tet.boud[2] = -1;
        tet.boud[3] = oldTetr.boud[3];
        vNewTetr.push_back(tet);

        tet.node[0] = m_vEdgeNode[vEdgeIndx[0]];
        tet.node[1] = m_vEdgeNode[vEdgeIndx[1]];
        tet.node[2] = m_vEdgeNode[vEdgeIndx[5]];
        tet.node[3] = m_vEdgeNode[vEdgeIndx[4]];
        tet.boud[0] = -1;
        tet.boud[1] = oldTetr.boud[1];
        tet.boud[2] = -1;
        tet.boud[3] = -1;
        vNewTetr.push_back(tet);

        tet.node[0] = m_vEdgeNode[vEdgeIndx[3]];
        tet.node[1] = m_vEdgeNode[vEdgeIndx[0]];
        tet.node[2] = m_vEdgeNode[vEdgeIndx[2]];
        tet.node[3] = m_vEdgeNode[vEdgeIndx[5]];
        tet.boud[0] = -1;
        tet.boud[1] = -1;
        tet.boud[2] = oldTetr.boud[2];
        tet.boud[3] = -1;
        vNewTetr.push_back(tet);

        tet.node[0] = m_vEdgeNode[vEdgeIndx[0]];
        tet.node[1] = m_vEdgeNode[vEdgeIndx[1]];
        tet.node[2] = m_vEdgeNode[vEdgeIndx[2]];
        tet.node[3] = m_vEdgeNode[vEdgeIndx[5]];
        tet.boud[0] = oldTetr.boud[0];
        tet.boud[1] = -1;
        tet.boud[2] = -1;
        tet.boud[3] = -1;
        vNewTetr.push_back(tet);
    }
}

// optimize the mesh
void RefineTetr::fn_OptimizTetMsh(std::vector<Node> &vR, std::vector<Tetr> &vElem)
{
    Node norm, e1, e2;
    Real dot;
    int m;

    for (auto &tet : vElem) {
        e1.x = vR[tet.node[1]].x - vR[tet.node[0]].x;
        e1.y = vR[tet.node[1]].y - vR[tet.node[0]].y;
        e1.z = vR[tet.node[1]].z - vR[tet.node[0]].z;
        e2.x = vR[tet.node[2]].x - vR[tet.node[0]].x;
        e2.y = vR[tet.node[2]].y - vR[tet.node[0]].y;
        e2.z = vR[tet.node[2]].z - vR[tet.node[0]].z;
        norm.x = e1.y * e2.z - e1.z * e2.y;
        norm.y = e1.z * e2.x - e1.x * e2.z;
        norm.z = e1.x * e2.y - e1.y * e2.x;
        dot = sqrt(norm.x * norm.x + norm.y * norm.y + norm.z * norm.z);
        e1.x = vR[tet.node[3]].x - vR[tet.node[0]].x;
        e1.y = vR[tet.node[3]].y - vR[tet.node[0]].y;
        e1.z = vR[tet.node[3]].z - vR[tet.node[0]].z;
        dot = (norm.x * e1.x + norm.y * e1.y + norm.z * e1.z) / dot;
        if (dot < 0.0) {
            m = tet.node[1];
            tet.node[1] = tet.node[2];
            tet.node[2] = m;

            m = tet.boud[1];
            tet.boud[1] = tet.boud[2];
            tet.boud[2] = m;
        }
    }
}
